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The description of the inspiral of a stellar-mass compact object into a massive black hole sitting 
at a galactic centre is a problem of major relevance for the future space-based gravitational-wave 
observatory LISA (Laser Interferometer Space Antenna), as the signals from these systems will be 
buried in the data stream and accurate gravitational-wave templates will be needed to extract them. 
The main difficulty in describing these systems lies in the estimation of the gravitational effects of 
the stellar-mass compact object on his own trajectory around the massive black hole, which can 
be modeled as the action of a local force, the self-force. In this paper, we present a new time- 
domain numerical method for the computation of the self-force in a simplified model consisting of a 
charged scalar particle orbiting a nonrotating black hole. We use a multi-domain framework in such 
a way that the particle is located at the interface between two domains so that the presence of the 
particle and its physical efi^ects appear only through appropriate boundary conditions. In this way 
we eliminate completely the presence of a small length scale associated with the need of resolving the 
particle. This technique also avoids the problems associated with the impact of a low differentiability 
of the solution in the accuracy of the numerical computations. The spatial discretization of the field 
equations is done by using the pseudospectral collocation method and the time evolution, based on 
the method of lines, uses a Runge-Kutta solver. We show how this special framework can provide 
very efficient and accurate computations in the time domain, which makes the technique amenable 
for the intensive computations required in the astrophysically-relevant scenarios for LISA. 
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I. INTRODUCTION 

One of the main sources of gravitational radiation for 
the future ESA-NASA gravitational wave observatory, 
the Laser Interferometer Space Antenna (LISA) [l], is 
the capture and inspiral of a stellar-mass compact ob- 
ject (SCO), with masses in the range m = 1 — 50Mq, 
into a massive black hole (MBH) in a galactic centre, 
with masses in the range AI = 10* — 1O^M0. These 
systems are usually known as Extreme-Mass-Ratio In- 
spirals (EMRIs) since the mass ratios involved are in the 
range /i — ra/M ^ 10^*" — 10~^. During the inspiral 
phase the system is driven by the emission of gravita- 
tional radiation, and hence there is loss of energy and 
angular momentum that makes the orbit shrink until the 
SCO plunges into the MBH. It is expected that LISA 
will be able to detect 10 - 10^ EMRI/yr |2] |3] up to 
distances within 1 Gpc [1] (see also for more details 
on the astrophysics of EMRIs). Recently, it has been 
suggested (6] [3 that inspirals of SCOs into intermediate- 
mass black holes, with masses in the range 50 — 350 Mq 
and presumably located in globular clusters, could be 
detected by future second-generation ground interferom- 
eters Hke Advanced LIGO |8] and Advanced VIRGO gj. 
The mass ratios are in the range 10"^ — 10"^, and in 
consequence they are called Intermediate-Mass- Ratio In- 
spirals (IMRIs). It is expected that techniques to de- 
scribe EMRIs may also been used for IMRIs, at least to 
a certain degree of precision. 

During the long inspiral, an EMRI will spend many 



cycles inside the LISA band, of the order of ~ 10^ during 
the last year before plunge into the MBH [10]. However, 
the gravitational-wave signals from EMRIs will be buried 
in the LISA data stream with the instrumental noise and 
the gravitational wave foreground (produced by compact 
binaries in the LISA band). To extract these signals and 
the relevant physical parameters that characterize them, 
we need to have a very precise a priori theoretical knowl- 
edge of the gravitational waveforms. This means to de- 
scribe the inspiral of the SCO taking into account the 
gravitational backreaction, that is, the influence of the 
SCO gravitational field on its own motion. This is an in- 
teresting but difficult theoretical problem, and different 
methods have been developed to solve it (see |111[T2][T3| ). 

While techniques for constructing templates good 
enough for detection are getting ready, mainly based on 
the use of the adiabatic approximation |14| [T5| [TBI 113 1 
methods to build templates good enough for extraction 
of physical information are not yet fully developed. The 
main difficulty being that one requires a more precise 
treatment of the self-gravity of the SCO and its impact on 
the gravitational waveform. In relation to this fact, there 
is currently a significant activity on the study of the ac- 
curacy of the different types of adiabatic approximations 
that have been introduced [HI [HI EOl EH El- On the 
other hand, the self-force approach to the equations of 
motion |23l [23] (see also jlll [25l [26] ) is an step forward 
to a precise estimation of the radiation-reaction effects 
and in consequence, towards the construction of accurate 
waveform templates. In this approach the backreaction 
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effects on the SCO are described as the action of a lo- 
cal force, the self-force, which can be computed in terms 
of the perturbations generated by the SCO with respect 
the MBH background spacetime. In practice, to compute 
the self-force one needs to regularize the perturbations, 
similarly as it happens in electromagnetism |27l [28] . To 
that end, a mode sum regularization scheme has been de- 
signed (in the gravitational case it has been formulated in 
the Lorenz gauge) |2l|30l|3l]|3l|33ll3lll35l|3g. It tells 
us how to subtract, mode by mode (for a Schwarzschild 
MBH the modes correspond to a harmonic decomposition 
of the perturbations), the singular part of the perturba- 
tions that does not contribute to the self-force. There- 
fore, what we need is a method of computing the full 
retarded solution of the perturbative equations for ap- 
plying the mode-sum regularization scheme and obtain 
in this way the self-force. The perturbative field equa- 
tions are a set of ten linear partial differential equa- 
tions (PDEs) for the metric perturbations h^^ {h^^ — 
g/xi/ ~ i/ii/ linear order, where g^^^ is the spacetime 
metric and g^^ is the background metric describing the 
MBH), which only in certain gauges (e.g., the Regge- 
Wheeler gauge |33) can be decoupled. In any case, to 
solve them completely we need to resort to numerical 
methods. Following the initial studies of black hole quasi- 
normal modes, frequency methods were used successfully 
[381 IMl iQl mi mils], and it was found they provide ac- 
curate results for EMRIs with moderate eccentricities. 
However, the frequency domain approach has more diffi- 
culties with highly eccentric orbits, which are of interest 
for LISA, since one has to sum over a large number of 
modes to obtain a good accuracy, and convergence may 
be an issue. This has opened the door to time-domain 
methods, which are not affected much by the eccentric- 
ity of the orbit and may be more efficient for the case 
of high-eccentricity EMRIs. In the last years there has 
been an intense activity on this front, both for a nonro- 
tating background |36l gH iH miZl gS] iS] , and for a 
rotating background |50l [5T] |52] . The main drawbacks of 
time-domain methods have mainly two origins: (i) The 
fact that one has to resolve very different physical scales 
(both spatial and temporal) present in the problem due 
to the extreme mass ratios involved (see, e.g. [53]). That 
is, using a standard numerical discretization of the prob- 
lem we are led to resolve the typical gravitational wave- 
lengths (comparable to the size of the MBH) and, at 
the same time, scales in the vicinity of the SCO, which 
are crucial for evaluating the self-force. This translates 
in a demanding requirement of computational resources, 
(ii) The fact that the SCO is described as a point-like 
object. This introduces Dirac delta distributions in the 
SCO energy-momentum distribution that lead to loss of 
differentiability in the solution of the perturbative field 
equations. This fact can degrade the convergence prop- 
erties of the numerical algorithms used. Moreover, such 
a localized distribution of matter can also introduce spu- 
rious high-frequency modes that contaminate the numer- 
ical solution and, in consequence, degrade its accuracy. 



Recently, there have been different proposals to im- 
prove the performance of time domain methods. Barack 
and Goldbourn |5l] have introduced a new technique to 
compute the scalar field generated by a pointlike scalar 
charge orbiting a black hole. This technique consists in 
subtracting from each azimuthal mode (in the Kerr ge- 
ometry the field equations are not fully separable in the 
time domain and one has to tackle them in 2+1 dimen- 
sions) of the retarded field a piece that describes the sin- 
gular behavior near the particle. This is done through 
a careful analytical study of the scalar field near the 
particle, using a puncture scheme which resembles the 
puncture model used for simulations in numerical rela- 
tivity |55l l56) [57] . This technique has been extended to 
the electromagnetic and gravitational cases by Barack, 
Goldbourn and Sago [58]. On the other hand, Vega and 
Detweiler jl^ have introduced another new method for 
regularizing the solution of the field equations. Their 
approach, tested on a simplified model of a charged par- 
ticle orbiting a nonrotating black hole, regularizes the 
retarded field itself by identifying and removing first, in 
an analytical way, the singular part of the retarded field. 
This alternative approach to the mode-sum regulariza- 
tion scheme yields a finite and different iable remainder 
from which the self-force can be computed. This remain- 
der is the solution to a field equation with a nonsingu- 
lar source, which avoid the problem (ii) above. Finally, 
Lousto and Nakano [59] have also introduced an analyt- 
ical technique to remove the particle singular behaviour. 
Their method is global and also produces a well behaved 
source for the field equations. 

Whereas these new techniques help in dealing with 
problem (ii) above, they do not completely solve the 
problem (i) since the regular source terms that these 
new schemes produce still have associated with them a 
length scale (or, from the numerical point of view, there 
are still special spatial resolution requirements associated 
with those source terms). In this paper, we introduce 
a new time-domain scheme towards the computation of 
the self-force which, for the case of a nonrotating black 
hole, eliminates completely any length scale associated 
with the SCO. This is done by using multiple subdo- 
mains and locating the particle in the interface between 
two of them (this has similarities with what was done 
in [16] using the finite element method, where in one of 
the numerical schemes proposed the particle was located 
between two elements). In this way, the Dirac delta dis- 
tributions do not appear in our equations and the pres- 
ence of the SCO enters through the boundary conditions 
that communicate the solutions at the different subdo- 
mains. As a consequence, we are solving wave-type equa- 
tions with smooth solutions, which avoid the problems 
described in (ii) (preliminary results have been reported 
in pCT). Regarding (i), we just need to provide the nu- 
merical resolution to describe the field near the particle, 
but not the particle itself, which makes the computation 
much more efficient. Our numerical algorithms are based 
on the PseudoSpectral Collocation (PSC) method (see. 
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e-g- IfiT]), which has been applied to numerical relativ- 
ity [62], and recently it has also been used in |63] for 
one-dimensional head-on collisions of black-holes. And 
very recently, in |64j . a discontinuous Galerkin method 
has been introduced and some gravitational waveforms 
for extreme-mass-ratio binaries are computed. This work 
uses similar techniques to the ones introduced in [igl and 
the ones that we present here. 

In this paper, we describe a set of techniques and meth- 
ods to use the PSC for the computation of the self-force 
on a charged scalar particle in circular orbits around a 
non-rotating black hole. We also show some results of 
the numerical implementation. The organization of the 
paper is as follows: In section |ll] we introduce the ba- 
sics of the model of a charged scalar particle orbiting a 
nonrotating black hole, including the basic formulae for 
the computation of the self-force via the mode-sum regu- 
larization scheme. In section ITTll we introduce all the in- 
gredients of a new time-domain numerical framework for 
the computation of the self-force in such scenario, from 
the mathematical foundations to the practical implemen- 
tation of the computations. In section |IV| we show the 
performance of a numerical code we have designed to im- 
plement the new scheme, and results of the computation 
of the self-force, in particular for the innermost stable 
circular orbit. In section |V] we draw conclusions from 
the results shown and discuss possible future avenues in 
the development of these techniques for the simulations 
of EMRIs in relevant physical situations. Throughout 
this paper we use the metric signature (— , -|-, -|-) and 
geometric units in which G = c = 1. 

II. DESCRIPTION OF THE MODEL: 
CHARGED SCALAR PARTICLE ORBITING A 
NONROTATING MBH 

In this paper we present a new technique for the com- 
putation of the self-force. In order to simplify things we 
focus in the particular case of a charged scalar particle 
in circular geodesies around a non-rotating MBH. In this 
simplified model, the spacetime metric is not dynamical, 
in the sense that it is not affected neither by the parti- 
cle nor by the scalar field that it generates. In our case, 
since we are considering non-rotating BHs, the metric is 
the Schwarzschild metric, which can be written as fol- 
lows: 

ds^ = f{-dt^ + dr*^) + r^dn^ , dn^ = d9^ + sin^ 0dip^ , 

(1) 

where (x^) ~ {t, r, 9, (p) are the so-called Schwarzschild 
coordinates, /(r) — l — 2M/r (where M is the BH mass), 
and r* is the tortoise coordinate, given by: 

r*^. + 2Afln(^-l). (2) 

In this geometry we assume there is a particle with scalar 
charge q, associated to a scalar field ^{x'^). This particle 
is orbiting the BH, and in doing so generates scalar field 



which in turn influences the particle trajectory. That 
is, the particle motion is affected by the field created by 
itself. In this way, this model contains all the ingredients 
of the gravitational case, in which the particle motion is 
influenced by its own gravitational field. 

The equation for the scalar field is then (see, e.g. jllj): 

g"''V„V^<i>(x) = -Airp = -Anq f dr S^{x, z{t)) , (3) 

that is, a wave-type equation with a source term that 
describes the particle energy density due to its scalar 
charge. In this equation, the spacetime metric g^^ is 
the Schwarzschild metric denotes the associated 

canonical connection; r denotes proper time associated 
with the particle along its timelike worldline 7, which we 
denote by x^ — z^{t); 64^{x,x') is the invariant Dirac 
functional in Schwarzschild spacetime, which is defined 
by the relation 

/ d^x./^^f{x)5,{x, x') = fix') , (4) 

and the equivalent one for the other argument. In this 
relation, g denotes the metric determinant. Taking into 
account the properties of 64^{x,x'), it follows that the 
source term in the scalar field equation ([3} only has sup- 
port on the particle worldline 7. 

The equations of motion for the particle trajectory 
{x'^ = z'^{t)) that one would obtain from energy- 
momentum conservation are: 

m—^F^^q{gf"' + u''u'')V,'i>, uf"^—, (5) 
dr dr 

where m and m'^ are the particle mass and 4-velocity, 
respectively. However, this expression is a formal one 
due to the fact that the force diverges on the particle 
worldline 7 (see [65] for a derivation of the regularized 
equations of motion). An analysis of the solutions of ([3} 
and ([5J reveals (see for details ^T\) that the gradient 
of the field, V^<i>, can be split into two parts [25]: A 
singular piece, <i>^, which contains the singular structure 
of the field and satisfies the same field equation, that is, 

g^^V^V^*^ =-4^p, (6) 

and a regular part, <i>^ = <i> — <i>^, which satisfies an ho- 
mogeneous wave equation 

g"'^v„v^a>- = , (7) 

and which is solely responsible of the deviation of the 
particle from geodesic motion around the BH, 

m— = g(g^''-Ku^unV^$'^. (8) 
dr 

This part is associated with the tail part of the scalar 
field. This is the analogous equation to the MiSaTaQuWa 
equation of the gravitational case (see [2F, i24j), and = 
q(gPi'_l_yPyi')y^(j)R is called the self-force on the particle. 
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In order to solve the equation for the scalar field 
[Eq. (|3j] it is very convenient to take advantage of the 
spherical symmetry of the Schwarzschild spacetime and 
decompose <I> in scalar spherical harmonics, Y^{6, (p) (see 
Appendix All, 



(9) 



where the harmonic numbers {£, m) take the usual val- 
ues: e = 0,1,2, ...,cx), and m = -£,-£+ 
Since we want to compute the self-force on the parti- 
cle for geodesies orbits, and since geodesic motion takes 
place on a plane in Schwarzschild geometry, we will as- 
sume, without loss of generality, that the plane is given 
hy 9 = tt/2. We will also parameterize the motion of 
the particle in terms of the coordinate time t, instead of 
proper time r. That is, the particle world-line, 7, will be 
given by {t, rp{t), Tr/2, (pp{t)) . Taking this into account, 
we can introduce the expansion ^ into the scalar field 
equation ([3} and find that the equations for the differ- 
ent harmonic coefficients ^™{t, r) decouple and have the 
form of a 1+1 wave-type equation: 



Or* 



where Vg{r) is the Regge- Wheeler potential for scalar 
fields on the Schwarzschild geometry, given by 



Veir) = fir) 



£{£+!) 2M 



(11) 



and S™ is the coefficient of the singular source term gen- 
erated by the particle: 



(12) 



where an overbar denotes complex conjugation. 

On a hypersurface {t — to} we can prescribe initial 
data for $7, {o^f.o^f) = (^^ (*o, r), 9t$7'(to, r)), and 
then find the corresponding solution (given that <i> satis- 
fies a wave equation, the problem is well-posed). There 
are a couple of comments in order: First, the solution 
found in this way corresponds to the (^, m) contribution 
to the full retarded solution. Second, this solution will be 
finite and continuous at the particle location, but it will 
not be different iable at the particle location in the sense 
that the radial derivative from the left and from the right 
of the particle yields different values. Moreover, the sum 
of the multipole coefficients over £ will diverge at the par- 
ticle location. This can be fixed by extracting, multipole 
by multipole, the singular part of the scalar field. Since 
we are interested on regularizing the self-force, which is 
defined at the particle location, let us introduce first a 
multipolar decomposition of the gradient of the scalar 
field 



E 



<^T{t,r)Yr{e,ip). 



(13) 



so that the gradient of the retarded field is given by 



1=0 



(14) 



Obviously, ^^{x^) also diverges at the particle world- 
line although the <i>^ are finite. Following the discussion 
of the previous section, the gradient of <I>, and hence the 
self-force, can be regularized by splitting the full retarded 
field into a singular part and a regular part (see, for in- 
stance, jlll 125) 1. such that they satisfy equations (|6|-(|8|. 
That is, the regular part of the gradient of the scalar field 
will be given, on the particle location, by 



00 



<i>i^(^(r))= lim 



{¥^{xn-^'/{^n) ■ (15) 



=0 



The singular field is known in a neighborhood of the 
particle worldline. In particular, the multipoles of the 
singular part of the gradient of the scalar field at the 
particle worldHne are given bv [29 l l30l ISTl l32l l33] l34l l66] : 



lim ^f''^ 

a:^ — >zi^ (r) 



(2£- l)(2f + 3) 



(16) 



where A^^ B^, C^, D^, ... are called the regulariza- 
tion parameters |23l |33] \M[ 166] . They are independent 
of £, but depend on the particle dynamics. The singu- 
lar part corresponds to the first three terms which lead 
to quadratic, linear, and logarithmical divergences when 
we sum over £. The remaining terms form a convergent 
series that does not contribute to the self-force (each of 
them). In our approximate calculations, we maintain the 
term as it accelerates the convergence of the series as 
we increase the number of multipoles included |35l 136] . 
For the case of interest of this work, circular geodesies, 
the non- vanishing regularization parameters are A^., B,., 
and j29t i34t l35] . which can be written as follows: 



A^ 



' rl l-2M/r 



B^ 



1 l-MI/r^ 



2M/rp 



1 - 3M/rp 
^1/2 - 2(l-2Af/rJ^3/2 



(17) 



(J8) 



D,. = 



\_ 2{l-2M/Tp) 
rl\ l-iM/Vp 
(l-M/rp)(l~4M/r^) 

8(l-2M/rp) "1/2 
(1 - 3Af /rp) (5 - 7M/rp - WM^/r^) 



M 1 



3M/rp "1/2 



16(1 - 2M/rpY 

3(l-3Af/rp)2(l + Af/rp) 
16(1 - 2M/r )2 



3/2 



5/2 



(19) 
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where Up is a sign that takes the value +1 when the limit and 



in equation (16 I is performed from the right (r > r^) and 



— 1 when thehmit is performed from the left (r < Tp). 
The coefficients Fq can be computed in terms of hyperge- 
ometric functions as follows (see [67], equation (15.1.1)): 



2^1 



Q,^;i; 



M 



rpf{rp) 



(20) 



Since the only non-vanishing component of the regular- 
ization coefficients is the radial one, this is the only com- 
ponent of the self-force that is actually singular, and 
hence the only one to be regularized. In this way the 
regularized self-force is: 



F^^q^^iz'^ir)). 



(21) 



III. 



A NEW TIME-DOMAIN FRAMEWORK 
FOR SIMULATIONS OF EMRIS 



In this section we introduce all the ingredients of a 
new time-domain method to evolve the equation(s) for 
a scalar field generated by a charged particle orbiting a 
nonrotating black hole, and also to compute the self-force 
from the evolved solution. 



A. Mathematical Preliminaries 

We are going to introduce some mathematical devel- 
opments needed in order to compute, using numerical 
algorithms based on the PSC method, the self-force on 
a charged particle on circular geodesies. We discuss first 
the particular formulation that we implement numeri- 
cally. Since our numerical framework deals with multiple 
domains it is useful to adopt a first-order formulation of 
the scalar field equation. The main reason is that first- 
order formulations of PDEs can be adapted to the hy- 
perbolic character of the underlying equation, and hence 
they are suitable for a correct communication between 
domains. 



We perform a first-order reduction of ( 10 I by introduc- 
ing the following new variables 



(22) 
(23) 
(24) 



The evolution equations for U = {ijj'^ , , (p™) , that is 
for a given harmonic (£, m), written in matrix form, are: 



where 



dtU 




A-d^-M + E 




-s, 


(25) 











1 0\ 
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0, 


(26) 




1 






0/ 





cm 

o,-^^(.*-.;w),o 



(27) 



It can be seen that (25 I is a first-order symmetric hyper- 
bolic system of PDEs, that is, it has the maximum degree 
of hyperbolicity, as expected of a system that comes from 
the reduction of a wave-type equation. The characteristic 
structure is described by the matrix A. 

Given that discontinuities in the solution are only al- 
lowed across characteristics, and given that due to the 
singular character of the source some of our variables 
will have jumps across the radial particle location, it is 
important to study these jumps using the formulation 
just introduced. To approach this question it is conve- 
nient to divide the spatial domain (the radial direction 
as parametrized by the coordinate r*) into two disjoint 
regions, one to the left of the particle (r* < r*{t)), and 
one to the right of the particle (r* > r*(t)). Then, we can 
write the solution of the equations ( [25| in the form (see 
also 06]): 

U = U_{t, r)e{r;{t) - r*) + U+{t,r)Q{r* - r;{t)) (28) 

where Q denotes the Heaviside step function. Introduc- 
ing (28 1 into (25 I we can derive the jump across the par- 



ticle location, defined as 

[Al - lim \^{t,r*)- lim^A_(t,r*), (29) 

P P 

for the different variables of our problem. We find the 
following result: 



0, [^T 



sv 



firp) 



(30) 



These are the conditions we have to impose in our nu- 
merical evolution algorithm. In practice, these condi- 
tions have to be imposed on the characteristic fields, 
the fields associated to the characteristics. In our prob- 
lem there are two different types of characteristics: (i) 
The {t — const.} surfaces, and (ii) the null surfaces 
{t ± r* — const.}. The associated characteristic fields 
are -0™, and 0™ =F (p™, respectively. 

Finally, our system of equations has to be comple- 
mented with initial conditions and boundary conditions. 
The initial conditions, since we are dealing with a first- 
order system of PDEs, consist in prescribing U^{r*) = 
U{t = t^,r*) . We need boundary conditions at the ends 
of the spatial domain. The physical spatial domain is 
r* G (— oo, -l-cx)), with r* — > — oo corresponds to the 
horizon location, while r* — > +oo corresponds to spa- 
tial infinity. In our numerical computations we will take 
a truncated domain: r* g [?'h7''i*]j ^ind therefore we need 
to prescribe outgoing boundary conditions at the ends 
of that domain (also called absorbing boundary condi- 
tions). As an approximation we take the conditions that 
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are exact for the case in which there is no potential, that 
is 



(31) 
(32) 



This is the leading approximation for the outgoing 
boundary conditions. They can be improved by analyz- 
ing the solution near r* ±oo. However, we can al- 
ways take values of (r*,r*) such that the boundaries are 
not in causal contact with the particle location, avoid- 
ing contamination of the solution due to propagation of 
unphysical modes from the boundaries. 



B. Using the PseudoSpectral Collocation Method 

We now describe in some detail the numerical tech- 
niques we use to solve for the full retarded scalar field 
and to compute from it the regularized self-force. 

In order to solve for the PDEs that describe the scalar 
field, equations (25 1, we use the PSC to discretize in 
space 



Once this is done we obtain a system of ordi- 
nary differential equations (ODEs) that can be solve by 
using the method of lines (see, e.g. |68]) applying a con- 
venient ODE solver. In general, spectral methods can 
approximate solutions of PDEs by expanding the vari- 
ables in a given basis of functions, and then by using an 
appropriate criterium that forces this expansion to ap- 
proach the exact solution as we increase the number of 
functions included in the expansion. In the case of the 
PSC method, the criterium consist in imposing the solu- 
tion exactly at a set of collocation points (see, e.g. |61j , 
for details on the PSC method). In our work, we use the 
Chebyshev polynomials, {T.^{X)} {X e [-1,1]), as the 
basis functions (see Appendix A 2 for definitions), and 



for the collocation points we take a Lobatto-Chebyshev 
grid. If we want to use a grid of N collocation points, 
the Lobatto-Chebyshev grid is made out of the zeros of 
the following polynomial (see Appendix A 2 I : 



(33) 



where the prime indicates differentiation with respect to 
X. The zeros can be written as follows: 



X, 



(i = 0,l,...,7V) 



(34) 



Now, we need to map the domain in which the Cheby- 
shev polynomials are defined (we will call it the spec- 
tral domain) to the physical domain. Before doing that, 
we have to mention that for computational reasons that 
we discuss later, we split the computational domain, 
— [r'^,r}], into a number of subdomains {D): 



(35) 



where r* ^ and r* ^ are the left and right boundaries of 
the subdomain U^. They are disjoint subdomains, that 
'^^ ^a-i R~'''aL (^^^ figure [T] and the caption there). 



Particle location 
(the nodes are identified) 



N N 



N 



H 



a,R 'a+1,L 



FIG. 1: The figure shows the structure of the one-dimensional 
spatial grid, the division in subdomains and the location of 
the particle at the interface between two of them. The bound- 
aries of the domain have coordinates (this coordinate will 
have a finite value and is meant to approach the BH horizon 
at r* = — oo) and Vj (this coordinate will have too a finite 
value and is meant to approach spatial infinity at r* = -l-oo). 
The particle is at the boundaries of two different boundaries 
(which are identified), which have coordinates r* ^ and r*^i j^ 
that satisfy the relation (471. 



We apply the PSC method to each subdomain, in 
the sense that our variables have different expansions in 
Chebyshev polynomials in each subdomain. These differ- 
ent expansions are related by using appropriate boundary 
conditions that we discuss in section [III C[ To apply the 
PSC method to each subdomain, we map the physical 
subdomain ^la (a = 1, ■■•,£') to the spectral domain, 
[-1,1], using a linear mapping X^ : [r*_i,r* ^ — > 
[-1, 1] , with 



2r* 



(36) 



' a,R 



The inverse correspondence is another linear mapping. 



[-1,1] — Kl.^a] - with 



X 



' a.R 



X 



' a,L 



' a,B. 



(37) 



The variables of our problem are arranged in a vector, 
U. Then, at a given domain {a — 1, ... , D), we have 
the following spectral expansion of our variables: 

N 

C/Ar(i,r-*)-^a„(t)r„(Xjr*)), (38) 

where the a„ are (time-dependent) vectors that contain 
the spectral coefficients of the expansion of our variables. 
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In the PSC method, we have also a physical expansion, 
which looks as follows: 



N 



(39) 



i=0 



where C^{X) are the cardinal functions |6T] associated 
with our choice of basis functions (Chebyshev polyno- 
mials) and set of collocation points (Lobatto-Chebyshev 
grid) [see Appendix A 2 for details]. The cardinal func- 
tions have the following remarkable property: 



(40) 



In this way, the time-dependent (vector) coefficients, 
{J7j}, of the expansion (39 1 are the values of our vari- 



ables at the collocation points 



U{t,r*{X,))^U,{t) 



(41) 



These are the variables that one looks for in the PSC 
method. The spectral (38 1 and physical (|39jl represen- 



tations are related via a matrix transformation (see Ap- 
pendix A 2 I. The computations (float-point operations) 



required to change representation scale, with the number 
of collocation points, as iV^, as it can be deduced from the 
fact that it is a matrix transformation. However, using 
the change of spectral coordinate given in equation ( A8 1, 



we can perform the change of representation by means of 
a discrete Fourier transform using a Fast-Fourier Trans- 
form (FFT) algorithm. In our numerical codes, we use 
the routines of the FFTW library j69]. Then, the number 
of computations required for a change of representations 
scales as TVlniV with the number of collocation points. 

Changing between representations is useful in order to 
perform some operations. For instance, differentiation 
is easier in the spectral representation, so we can trans- 
form from the physical to the spectral representation, 
compute derivatives there, and finally transform back to 
the physical representation. In the case of a Chebyshev 
PSC method, the differentiation process can be described 
by the following scheme 



9.. 



1-^, . FFT f . 9, 



{bj ^ {(a,.J7)J,(42) 



where are the spectral coefficients associated with 

the spatial derivative, and their relation to the coeffi- 
cients of the variables, {a„}, is given by (see, e.g. |61) ) 



bN = bN 



1 = 0, 



-{2na„ + b„+i} {n = N-l, 



and 



2 for n = , 
1 otherwise . 



(43) 
.1)(44) 

(45) 



In the PSC method, we find a discretization of our 



system of equations ( 25 I by imposing them at every col- 
location point. In practice, this is done by introduc- 



ing the expansion (39 1 into the equations (25 I, and then 



we evaluate the result at every collocation point of our 
Chebyshev-Lobatto grid ( [34| . We obtain a system of 



ODEs for the variables 



U, = A-id,,U),+M-U, + S, 



(46) 



where the dot denotes differentiation with respect to the 
time coordinate t, and {d^,U)^ has to be interpreted ac- 
cording to the scheme in equation (42 I. 



C. Evolution Algorithm 

In this section we discuss the details of the time evolu- 
tion of the discrete equations we derived in the previous 
section. In particular, we describe how to introduce the 
particle in the multi-domain PSC method that we pro- 
pose. The main argument is the following: If we put the 
particle inside one of the subdomains fi^ we are intro- 
ducing in the equations (251 a singular term that, in the 



case of a scalar charged particle, will produce solutions 
that are not differentiable (in the sense that the deriva- 
tive with respect to r* is not single- valued at the particle 
location). That is, our solution will not be smooth and 
hence we cannot expect the PSC method to converge ex- 
ponentially to the true solution. That would spoil one of 
the main motivations for using this numerical technique, 
that is, the accuracy that the PSC method provides. To 
avoid this (and this is part of the reason for using a multi- 
domain scheme) we put the particle at the interface be- 
tween two subdomains. If the subdomains are and 



f^a+iJ have: 



' aM 



-1,L 



(47) 



Since we have restricted ourselves to the case of circu- 
lar orbits, once the grid has been set we do not need to 
change it during the evolution. In the case of generic 
orbits, in order to maintain the particle at the interface 
between two subdomains we have either to make a co- 
ordinate change or implement a moving grid scheme (a 
similar situation was confronted in (1^). We discuss this 
further in section \V\ 

For each subdomain fi^ (a = 1, • ■ • , -D), we evolve the 
set of ODEs (46 1 independently. Since we are locating 



the particle at the interface between two subdomains, 
the last term in (46 1, S, does not appear. This term 



is the source term that accounts for the particle's en- 
ergy density [see equation (12l]. The main implication 
of this setup is that we have to solve the homogeneous 
field equations, which are equations with smooth solu- 
tions and hence, the advantages of the PSC method are 
preserved. Then, in our framework, the contributions of 
the particle to the solution appear as boundary condi- 
tions between the subdomains. In summary, we evolve 
the (homogeneous) field equations for each subdomain 
independently and connect their solutions by boundary 
conditions at the interfaces. The equations are evolved 
using a Runge-Kutta (RK) solver (see |701 FTT] '). typically 
a RK4 algorithm. 
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The key point in the evolution is the imposition of 
boundary conditions at the interfaces between subdo- 
mains. There are two possible situations for a given in- 
terface: (i) The particle is not there. In this case we only 
have to impose the continuity of the solution [that would 
correspond to imposing the junction conditions given in 
equation (30 1 with zero right-hand sides], (ii) The par- 
ticle is there. In this case we have to impose the junc- 
tion conditions in equation (30 1 where 5*™ is given in 



equation (12 1. The question is then how to impose these 



boundary conditions numerically within our framework. 
There are different ways of doing this, and in this work 
we have chosen to impose this boundary conditions at 
the interfaces between subdomains in a dynamical way, 
via penalty terms. 

The penalty method is a well-known technique and has 
been applied to several numerical schemes to solve PDEs 
(for the PSC method, see [72] and references therein). It 
is relatively simple to implement it for elliptic-type prob- 
lems, that is, for problems that do not require evolution 
in time. For time-dependent problems, it can also be im- 
plemented but the result is not always numerically sta- 
ble In our framework, the implementation of the penalty 
method goes as follows: First, let us consider the two do- 
mains around the particle's location, say il^ and ri^+ii 



that equation (47 1 holds. Then, let us consider the solu- 
tions at these two subdomains, U^^ and U^_^_l. The equa- 
tions for the inner points are just equations ([46} (with 
= 0), and the equations for the nodes at the inter- 
face between these subdomains, r* ^ and r*_|_j^ ^ (which 
are identified), are modified in the following way: For 
the subdomain fi^ (we have simplified the notation by 
dropping the harmonic indices i and to) 



a,R 



(48) 



Tj. 



K,R + Va.R 



dtVa,R = d^-'l^a^R- ^ {4'a,R + Va,R 

- {<t>a+l,L+^a+lx)-W]p} . (50) 

and for the subdomain ^a+i 

Ma+l 



a+l,L 



L - Va+l,L- '4, [V'a+1,L - V-a.fl] i (51) 

a+l,L 



-l.L 



fa+l,L - (<PaM 



^a,R)-W]p}, (52) 



dtVa+l,L = 9r,<l>a + l.L - {'/'a+l,L " "Pa+l.L 

- i^a,R-^a,R)-[v]p} , (53) 

where Vp = V{rp), the quantity [(pj^ is the analytic jump 



imposed by the dynamics and given by equation ( 30 1 , and 



V'a,i?,(i) = V'(i, , V'a+l,L(i) = '/'(i, rl^i,L) , (54) 
0a,fl(i) - Ht, tIr) , <Pa+iAi) = Ht, , (55) 

fa,R{t) = fit, > fa+l,L{t) = fi*, ''I+l.i) > (56) 



rameters. 



tions (48 1 



''i,',4:,v ^^'^ '^4>'!~<p,'v (constant) penalty pa- 
The structure of the penalty terms in equa- 



531 obeys the following rationale: The main 



idea behind of the penalty terms is to drive the dynam- 
ical system to satisfy a set of conditions that are not 
part of the original evolution equations (like constraints 
on the variables that have to be satisfied for all times 
or boundary conditions), and the strength of the driving 
(penalty) terms is controlled by the pen alty pa rameters 



a,R 



T,'V and t'^'^}''". In our case, in section III A 
tioned the fact that the junction conditions (301 have 
to be imposed on the characteristic field of our system 
of PDEs, and these characteristic fields and their associ- 
ated characteristic surfaces where given there. Therefore, 
we have constructed the penalty terms so that the junc- 
tion conditions that are satisfied are those corresponding 
to the characteristic fields, that is [restoring the to) 
indices]: 



Or i m I m 



± 



cm 

I 



(57) 



On the other hand, i(j can be seen as a subsidiary variable, 
in the sense that we can first evolve the equations for 
and ip and then, use the result to evolve ij;. Then, we 
can use a different way of imposing the continuity of tp, 
instead of the penalty method we can just replace the 
right-hand sides of equations (48 1 and (51 1 by 



dti^a+l,L = 



K,R 



*a+l,L)/2, 
*a+l.L)/2, 



(58) 
(59) 



which ensures the continuity of ip by construction and we 
have seen in our numerical experiments that it is numer- 



ically stable (see Sec. IV I. 

Regarding the global boundary conditions {near the 
horizon, r* — r*, and near spatial infinity, r* — r*), we 
impose them directly at the corresponding nodes, with- 
out using the penalty method (which is another option). 
Since we have made a first-order reduction of our origi- 
nal PDE (10 1, we have to adapt the outgoing/absorbing 
boundary conditions to the set of variables U = {4>, 4>, <f)- 
Then, the boundary conditions at r* — r* are 



Mix = -5tVl,L, 
Mix = d^'^L, 
and the boundary conditions at r* — r* are 



Mdm 
Mdm 
Mdm 



VD..R: 
M^t^D.R, 



(60) 

(61) 
(62) 

(63) 
(64) 
(65) 
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and we remark that the first subdomain is number 1 and 
the last one is number _D, the total number of subdo- 
mains. 



IV. RESULTS FROM THE SIMULATIONS 

In this section we describe the results from a series 
of numerical experiments that show the potential of the 
methods just proposed. To that end, we have developed 
a numerical code that implements the techniques previ- 
ously described. This code is based on the C language 
and the use of the GNU Scientific Library [73], mainly 
for calculations with special functions, and the FFTW 
library j69] for performing FFTs. In the implementation 
of the physical domain in the numerical code we have 
adopted a comoving tortoise coordinate: r* = r* —r*. In 
this way, the particle location (and we emphasize again 
that we only deal with circular orbits) is always r* = 0. 

The first test we have performed is to study the evo- 
lution of a simple wave equation using the formulation 
and methods described in the previous sections. This 
case corresponds to: £ ^ m = 0, M ^ (which implies 
= 0), and q = (that is, no source, 5™ = 0). The 
test consists in following the propagation of an initial 
Gaussian packet in a multidomain grid and to study the 
convergence of the numerical scheme as the number of 
collocation points per subdomain, N, increases. In Fig- 
ure |2] we show a convergence plot for a simulation that 
uses two subdomains, r* - r* e [-550 M, 0] U [0, 550 M], 
connected by the penalty method as described in sub- 
section |III C[ The truncation error has been computed 
in the subdomain where the Gaussian packet is present 
at a chosen time. As one can see, the truncation error, 
estimated as the absolute value of the last spectral co- 
efficient, |a^|, decreases exponentially with the number 
of collocation points, as expected in the PSC method for 
smooth solutions. 

The same convergence follows for an initial Gaussian 
wave packet propagating on a Schwarzschild background, 
that is, for the case in which only 5™ = 0. This is ex- 
pected as the mathematical structure of the equation is 
essentially the same. The additional ingredient is the ex- 
citation of quasinormal modes of the black hole by the 
initial wave packet. 

When we introduce the SCO, i.e. the particle, the 
situation is conceptually different. If we think in global 
terms, the presence of the particle implies that the global 
solution (the solution in the whole computational do- 
main) will not be smooth. Hence, we cannot expect expo- 
nential convergence for the solution of our problem. How- 
ever, as we have argued above, our multidomain frame- 
work avoids the presence of Dirac delta distributions in 
the equations by locating the particle in the interface be- 
tween subdomains. Then, the presence of the particle 
enters through boundary conditions, actually matching 
conditions between subdomains. Therefore, at each sub- 
domain we will have a smooth solution, and hence we 
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FIG. 2: This figure shows the dependence of the truncation 
error (estimated by the logarithm of the absolute value of the 
last spectral coefficient, logj^Q \aj^\, associated with the field 
tp) with respect to the number of collocation points. This 
error corresponds to an snapshot of the evolution of the clas- 
sical wave equation, for an initial condition given by a moving 
Gaussian packet. The data indicates the exponential conver- 
gence of the numerical method. 

expect our numerical solution to converge exponentially 
towards it. In our numerical experiments with a particle 
(remember we restrict ourselves to the case of circular 
orbits) , we take zero initial data, that is 

2PTito,r*) = ^Tito,r*) = vTito,r*) = 0. (66) 

This initial data is obviously not consistent with Ein- 
stein's equations, and as a consequence the evolution 
produces an initial unphysical burst. We have to wait 
until this unphysical feature propagates away in order 
to analyze the solution and obtain physically relevant 
results. In Figure |3] we show snapshots of the evolu- 
tion of a scalar charged particle in circular motion at 
the Last Stable Orbit (LSO), in principle the most de- 
manding case (in Schwarzschild the LSO is located at 
rp — 6 M) . The figure includes details of the different 
variables, {ipj^ , (j>™ , (p™) for £ = m = 2, near the parti- 
cle location. These snapshots illustrate the ability of our 
method to capture the structure of the solution near the 
particle, in particular the ability of resolving the jump in 
the radial derivative of the field (snapshot on the right in 
Figure |3j . 

In order to further validate our numerical code, we 
have performed simulations with the particle at the LSO 
changing the number of collocation points, while leaving 
the number of subdomains fixed. In this way we have 
checked that our method can achieve the exponential 
convergence in each individual subdomain. In Figure |4] 
we show a convergence plot obtained from simulations 
that use four subdomains: r* - r* e [-550M, -20M] U 
[-20 M, 0] U [0, 20 M] U [20 M, 550 M]. The number of col- 
location points, N, is the same at each subdomain, and 
in this way we have more resolution near the particle. 
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FIG. 3: We show snapshots of the evolution of the scalar charged particle in circular motion at the LSO for the mode I = m = 2. 
These simulations used 12 subdomains and 50 collocation points per subdomain. They show the evolution of the variables 
■i/>™ (left), (pj^ (center), and ip^ (right), after a substantial time has passed and a number of wave cycles have been generated. 
In particular they show how the jump in the radial derivative of the field $^ is resolved (this information is encoded in the 
variable </5™) in this multidomain computational framework. 



where it is most needed. The figure shows that indeed 
our numerical scheme has exponential convergence. 
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FIG. 4: This figure shows the dependence of the truncation 
error (logj^Q |ajv|) with respect to the number of collocation 
points for simulations of circular motion at the LSO (r^ = 
6M) corresponding to the I = m = 2 mode, in particular 
for the field ip2- The data, taken from the subdomain r* — 
G [0, 20 M] indicates the exponential convergence of the 
numerical method. 

The multidomain feature of our method is useful for 
another important reason, namely computational cost. It 
is not the same having N points in D domains that hav- 
ing N ■ D points in one single domain. This is an impor- 
tant fact to be considered in a situation where the need 
for resolution comes only from some isolated regions of 
the computational domain. This is what happens in our 
problem and the type of multidomain structure can have 
dramatic consequences in the calculations. First of all, 
for computations in a given time step, the first option in- 
volves less calculations. Second, from the evolution point 
of view, in the PSC method, the Courant-Friedrichs-Lax 
(CFL) condition on the allowed size of the time step, 
Ai, is more stringent than in Finite Differences (FD) 
schemes. For the PSC method the maximum allowed 
time step is AtcFL ""^ki? — ''"ll/i^^'^) (this is set by 
the minimum distance between two collocation points. 



which occurs at the boundaries of the subdomains [61]), 
that is, it goes like with respect to the number of 

collocation points, whereas in FD schemes it goes like 
1 /N. Then, dividing the domain in subdomains can help 
in having a bigger AtcpL- In addition, the multidomain 
scheme can be seen as a way of adaptivity, in the sense 
that we can construct small subdomains for the regions 
that need to be well resolved, essentially near the par- 
ticle, and large subdomains for the regions that do not 
need to have high resolution, essentially far away from 
the particle. 

In order to illustrate how the multidomain feature of 
our numerical framework works, we have performed sim- 
ulations with a fixed number of collocation points but 
changing the number of subdomains. The aim is to show 
how the solution improves by adding new subdomains. In 
our simulations, the most demanding region for resolu- 
tion is clearly near the particle. Apart from this we have 
the waves leaving the particle location and going both 
towards the horizon and towards spatial infinity. These 
waves are in principle easy to resolve, but for modes with 
high £ and m we have that the source term oscillates like 

exp{imfl.p{t—to)} (where fip ^ y^M/r^ is the coordinate 

angular velocity of the particle), and hence the wave- 
length gets reduced so that we have many more waves 
that for low m. These waves are moving away and need 
to be resolved. In Figure [5] we show snapshots of evolu- 
tions with 50 collocation points per subdomain, but with 
different number of subdomains, namely from 2 to 16 sub- 
domains (half of them to the left of particle and the other 
half to the right). The figure shows the mode i = 10 and 
m = 6 (so that the period of the waves is 7r/(3rip)) for 
the three fields {ip, cj), (p). We can see how for few subdo- 
mains the waves are unresolved, and that as we increase 
the number of subdomains the solution converges. There- 
fore, the multidomain structure is a very useful tool to 
achieve accurate results with a reasonable computational 
cost. The key point is to realize what is the optimal num- 
ber of subdomains, their size, and their distribution over 
the whole computational domain. Given that in our case 
the period of the waves changes with the harmonic num- 



11 



ber m, the optimal strategy for setting the subdomain 
will depend on it. However, since the main aim of this 
work is to show the principal ingredients of the method 
and to illustrate its performance, we have not explored 
the possibility of changing the subdomain structure as 
m changes. However, this is something that should be 
done for optimizing the computational cost in the case 
of systematic calculations of the self-force. In this sense, 
it would be convenient to have a deep understanding of 
the computational parameter space in order to automat- 
ically adapt these parameters to the physical parameters 
and optimize in this way the calculations. This would 
require an extensive investigation that will be done and 
presented elsewhere. 

The next step in the validation of the code is to 
compute the components of the self-force acting on 
the charged scalar particle moving in circular geodesies 
around a nonrotating black hole. This is the goal of this 
numerical scheme, that is, to provide accurate compu- 
tations of the self-force with reasonable computational 
cost. The calculation of the self-force require to compute 
the evolution of the real and imaginary parts [note that 
the source term in the evolution equation ( [To| is com- 
plex and hence both the real and imaginary parts are 
needed] of the {£,m) harmonic components of <i>™. Since 
<I> is a real scalar we have that for each {£, in) the rela- 
tion $7™ = (—1)™$™ holds and hence we do not need 
to compute the modes with m < 0. Given that we need 
to truncate the sum over {£, m) at a given £„ax, the total 
number of evolutions that we need to run for a single self- 
force calculation is iVe,oK,ti„„, = -I- -I- 2)/2], 
where here the brackets denote the the nearest integer 
to the argument. For the particular case — 20, 
a typical case in our calculations, we need to perform 
-^evolutions — 231 evolutlous. 

In Table [l] we present the results of the computations 
of the self-force vector (actually of the gradient of the 
regular scalar field) acting on a scalar particle in circular 
(geodesic) motion around a black hole at the following 
radii: r/M ^ 6, 7, 8, 10, 14, and 20. For these simu- 
lations we have used our multidomain framework with 
subdomains that contain 50 collocation points. We have 
used the same number of subdomains to the right and 
to the left of the particle. The size of these subdomains 
is typically Ar* — 20, specially those near the particle. 
As we get far away enough from the particle towards the 
boundaries (typically at r* — ±(500 — 700)M), the size of 
the subdomains is increased since we need less resolution 
in those regions. The total number of subdomains that 
we have used ranges from 12 — 34. We have observed 
that at some point increasing the number of subdomains 
(maintaining the size of r* = 20M) does not change sig- 
nificantly the results. Another important point to men- 
tion is the fact that we are computing the self-force right 
at the particle location. But since the particle is at the 
interface between two subdomains, we obtain two val- 
ues of the self-force components, one from the subdo- 



main to the left [let us call it fla as in equation (47l] 



^q' = ^aC'^a i?,)' ^^'^ the other one from the subdomain 
to the right (f^a+i), ^S'"^ = '^IK+i.l)- This also pro- 
vides us with a test of the numerical calculations as both 
values have to agree to a good gegree of precision. In 
Table |l] we show our results and compare them with two 
types of calculations in the literature: (i) Calculations 
based on a time-domain method that uses a characteris- 
tic formulation of the scalar field equations [36j; (ii) cal- 
culations based on a frequency-domain method [71]. As 
we can see, we can get a good numerical approximation 
to the self-force components using a modest amount of 
computational resources. It is important to mention that 
these computations have a relatively low computational 
cost (relative to the computational cost of time-domain 
simulations of this sort). The average time for a full self- 
force calculation of the type just described in a computer 
with two Quad-Core Intel Xeon processors at 2.8 GHz is 
always in the range 20 — 30 minutes. 



V. CONCLUSIONS AND DISCUSSION 

In this paper we have introduced a new time-domain 
technique for the simulations of EMRIs. The main ingre- 
dient of the method is to use a multi-domain framework 
in which the SCO, described as a point-like object, is lo- 
cated at the interface between two subdomains. In this 
way we have shown that this technique enjoys the ex- 
ponential convergence property of the PSC method and 
its accuracy and, at the same time, it is also an efficient 
method to make time-domain computations of the self- 
force. We have shown that we can achieve a good accu- 
racy in this computations by using a relatively low num- 
ber of collocation points. In this sense, the multidomain 
framework allows us to locate more collocation points in 
the region where they are most needed, i.e. around the 
particle, by choosing appropriately the number of sub- 
domains, and their size and location. Another positive 
property of our numerical scheme is that it can be eas- 
ily parallelized to be used in supercomputers, as we can 
assign the work of each subdomain to different CPUs. 
The only information that we need to communicate is 
the one necessary to satisfy the matching conditions be- 
tween subdomains. 

The calculations performed for this paper show the 
potential of this technique for the description of EMRIs. 
These results still leave room for improvement as we have 
not explore yet the full parameter space of the numerical 
method (number of collocation points per subdomain, 
number and distribution of subdomains, parameters of 
the spectral filter, penalty parameters, etc), which is wide 
enough. In this regard, it would be convenient to be able 
to adjust the computational parameters to the physical 
parameters of each mode in an automatic way. In addi- 
tion to this, there are several additions/modifications to 
the present numerical framework that may improve the 
present accuracy and efficiency of the calculations, some 
of which we will investigated in future calculations by the 
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FIG. 5: We show here three plots of an snapshot of the evolution of the scalar charged particle in circular orbital motion at the 
LSO for the mode ^ = 10 and m = 6. Each plot shows different evolutions of the variables tp (left), (j) (center), and ip (right), 
for a fixed number of collocation points (A'^ = 50), but for different numbers of subdomains {D = 2, 6, 10, 16). In this we can 
see how increasing the number of subdomains leads to a better accuracy (we can see how the waves get better resolved and the 
solution converges as we increase the number of subdomains) with a much less computational cost as if we would have increase 
the number of collocation points in a single domain computational framework. 



TABLE I: In this table we show the values of the regular field at the particle location, ($a'~, ^a'"^)", computed at ipp = 0. 
These results are obtained for for £„ax ~ 20, A'^ = 50, and D = 12 — 36. The minimum size of the subdomains, as measured 
in terms of the tortoise coordinate is Ar* = 20M, which corresponds to the subdomains near the particle location. The table 
shows our numerical results for different circular orbits (r/M = 6, 7, 8, 10, 14, and 20) and the relative errors with respect to 
the results obtained with a time-domain method in |47] . and with a frequency-domain method in |74j and |36j . 



r{M) 


Component 
of $S 


Estimation using 
the PSC Method 


Estimation from 
Frequency-domain 
calculations in |74l and |36J 


Error relative to 
Frequency-domain 
calculations in [74 1 and j36j 


Error relative to 
Time-domain 
calculations in 1471 


6 


(-I-?'-, -!>?■+) 


(3.60777,3.60778) • 10^"* 
(1.67364,1.67362) • 10"'' 
(-5.30422,-5.30438) • 10"=* 


3.609072 • 10-"* 
1.67728 • 10"* 
-5.304231 • 10"^ 


(0.03,0.03) % 
(0.2,0.2) % 

(4- 10-*, 10-3) % 


(0.12,0.12) % 
(0.18,0.18) % 
(6- 10-*, 10"=*) % 


7 


(-l'?'-,-l>^'+) 


(1.76638,1.76639) • lO""* 
(7.84007,7.84001) • lO"'^ 
(-3.2730,-3.2733) ■ 10"^ 


7.85067- 10"^ 


(0.13,0.13) % 




8 


(-l'f'-,$f'+) 
($?'-, $^'+) 


(9.76454,9.76457) • 10"^ 

(4.0835,4.0832) ■ 10^^ 
(-2.2115,-2.2108) • 10"^ 


4.08250 ■ lO"'' 


(0.02,0.01) % 




10 


(■fr'-.i.r) 

(-!-?■-,$?'+) 


(3.74362,3.74363) • lO"'^ 

(1.3804,1.3801) ■ 10"^ 
(-1.1860,-1.1858) • 10"^ 


1.37844 ■ 10"'' 


(0.14,0.12) % 




14 


(■fr'-.K^) 
($?■-, $^'+) 


(9.176912,9.176914) ■ 10"** 

(2.7269,2.7266) • 10"^ 
(-4.8383,-4.8389) ■ 10"* 


2.72008 ■ lO"'' 


(0.25,0.24) % 




20 


(-j>f'-,$f'+) 

($?'-, -!>?■+) 
($5'-,<E>^'+) 


(2.08859,2.08858) • lO"'' 
(4.9444,4.9440) ■ 10"^ 
(-1.92449, -1.92436) • lO""* 


4.93790 • 10"^ 


(0.13,0.12) % 





"We show the values of the gradient of the regular field instead of the components of the self-force for the sake of comparing with other 
results in the literature. 



present authors. We list here some of them: (i) Compact- 
ification of the computational domain. We can change 
the mapping between the spectral and physical domains 
incorporating the two infinities (the horizon location at 
r* —oo and spatial infinity at r* oo) into the cal- 
culation. In practice, this mapping only would need to 
be performed in the first and last subdomains. By do- 



ing this, the outgoing boundary conditions (32 1 would 



be exact instead of approximate as they are know. An 
alternative to this could be to improve the boundary con- 
ditions by analyzing the solution near the two infinities, 
like for instance in [75j in a frequency domain framework. 



(ii) Reduction of the time step. As we have mentioned be- 
fore, the PSC method has a more strict CFL condition 
as a FD scheme, namely At ~ N~'^ versus At ^ N~^. 
This is because of high density of collocation points near 
the boundaries and the fact that the minimum time step 
come from the minimum separation between points. A 
way of changing this is again to change the mapping be- 
tween the spectral and physical domains. One known way 
of doing this is the modification introduced by KoslofF 
and Tal-Ezer j7^, which was shown to lead to a time step 
restriction like in the FD case, that is At ~ N~^. (iii) 
Richardson extrapolation. In our calculations we have 
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computed harmonic components of the scalar field, <I>™, 
up to a certain ij^^x^ this work i^^ax = 15 — 20. In 
other words, we have truncated the expansion in spheri- 
cal harmonics. We can try to improve the final numerical 
results for the scalar field and derived quantities, like the 
self-force, by profiting from our analytical knowledge of 
the expansions in inverse powers of the harmonic number 
£. This can be done using numerical techniques based on 
the Richardson extrapolation, like in |35) . which can pro- 
vide an estimation of the different quantities of interest 

Beyond the case studied in this paper, circular orbits of 
a charged scalar particle, we will study in the future how 
to extend our time-domain computational framework in 
order to include generic (eccentric) orbits. The main dif- 
ficulty is obviously to deal with a particle moving the 
radial direction as our present framework assumes that 
it is located at a fixed value (circular motion). There are 
several ways in which we can attack this problem. One 
is to try to implement some type of moving grid tech- 
nique (as it was done in a similar situation in |l6]), but 
the reconstruction of the grid at every time step could 
increase significantly the computational cost of the evo- 
lution. Another possibility is to try to make a change 
of coordinates to coordinates comoving with the parti- 
cle, so that the techniques presented in this paper can be 
implemented in a straightforward way. 

On the other hand, we expect that MBHs sitting at 
galactic centers will have considerable spins {a/M ~ 0.7 
or bigger, where a is the Kerr spin parameter) and hence 
the MBH should be described by the Kerr metric instead 
of the Schwarzschild metric. This means to have a less 
symmetric background, instead of the spherical symme- 
try of Schwarzschild just the axisymmetry of Kerr. In 
that case one can separate the dependence on the az- 
imuthal angle and is left with 2-1-1 wave-type equations 
with singular terms. The main difficulty in this case is 
that each m-mode (coming from the separation of the 
azimuthal angle) diverges logarithmically at the particle 
location. Then, before trying to transfer the techniques 
presented in this paper to the Kerr case, one has to ap- 
ply before a regularization procedure as it has been done 
in |48l IMl |59] . Apart from this, we expect that most 
of the methods presented in this paper will be helpful in 
achieving efficient simulations of EMRIs in the case of a 
spinning MBH. 
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APPENDIX A: SPECIAL FUNCTIONS USED IN 
THIS WORK 

Here, we summarize the main conventions used for the 
special functions involved in calculations of this paper. 



1. Spherical Harmonics 

The expression we use for the scalar spherical harmon- 
ics is: 



where P™ are the associated Legendre polynomials [we 
use the same expressions as in ^7\, equations (8.6.6) and 
(8.6.18)] 

where ^ is a non-negative integer and m is an integer re- 
stricted to the following range: m G {~i , — 1 , . . . , £ — 
1,£). 



2. Chebyshev Polynomials 

The Chebyshev polynomials can be expresses as fol- 
lows: 



cos n cos 



(A3) 



and are defined in the interval [—1, 1] with |T„(Ar)| < 1, 
where n is the degree of the polynomial. Chebyshev poly- 
nomials are orthogonal in the continuum in the following 



(T T 



dX 



T^{X)T^{X) 



7rc„ 



(A4) 

where the coefficients c„ are given in equation (45 1. 

The set of collocation points for the spatial discretiza- 
tion of our PDEs by means of the PSC method are the 
extrema of the Chebyshev polynomials, together with the 
end points X — ±1 [that is, the zeros of the polynomial 
given in equation 



-1,1], a Chebys^ 



34 1]. These points form, in the interval 



lev-Lobatto collocation grid. The car- 



dinal functions associated with them are (i = 0, . . . , N): 



{l~X^)T^{X) 



{l-Xf){X-X^)T^{X 



(A5) 
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Once this set of collocation points is adopted, and tak- 
ing into account the properties of the Gauss-Lobatto- 
Chebyshev quadratures (see, e.g. [61]), the Chebyshev 
polynomials have another orthogonality relation, this 
time in the discrete, in the following sense {n, m = 
0,...,iV): 



N 



(A6) 



where are the weights associated with the Chebyshev- 
Lobatto grid, = tt/{Nc^), and where the q's are nor- 
malization coefficients given by 



2 ioT i^O,N , 
1 otherwise . 



(A7) 



Finally, the constants f„ in ( A6 1 are given by i/^^ — 7rc„/2. 

On the other hand, introducing a new variable, X = 
cos 9, the Chebyshev polynomials look like 



r„(cos6l) = cos(n6') . 



(A8) 



Then, an expansion in Chebyshev polynomials can be 
mapped to a cosine expansion. 



filter of the exponential type to the solution after every 
time step, that is, after every full RK step. The scheme 
for the action of the spectral filter is 



FFT 



Filter 



FFT 



{U.}, (Bl) 



where {iJj} are the values of the solutions at the colloca- 
tion points after the RK step; {a„} are their correspond- 
ing spectral components; are the filtered spectral 
components; and {C/j} are the filtered values of the so- 
lution at the collocation points. The exponential filter is 
defined by its action on the spectral coefficients {a„} to 
yield the spectral coefficients {b„}. This action is given 
by 



(B2) 



where a{n/N) is the exponential filter 



exp 



1 for < 71 ^ A?^, , 



(B3) 



APPENDIX B: SPECTRAL FILTER 

In order to reduce the spurious high-frequency com- 
ponents of our numerical solutions, we apply a spectral 



where is the cut-off mode number, 7 is the order of the 
filter (typically chosen to be of the order of the number 
of collocation points, N), and a is the machine accuracy 
parameter, which is related to the machine accuracy, e^^, 
by a = -Incjvf. 
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